# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating final dataset for Chicago.

#### libraries ####

library(MASS)
library(sp)
library(stargazer)
library(sf)
library(matrixStats)
library(tidyverse)
library(scales)
library(tidycensus)
library(spdep)
library(gridExtra)
library(readxl)

#### downloading census data #### 


# Download Decennial Census data 2010

census_api_key(key = 'f1ad9e98451897912bfdf1227c9a357eaa7d0d78', install = TRUE)

vars_census <- c("P001001", "P012006", "P012007", "P012008", "P012009", "P012010",
                 "P012011", "P012012", "P005003", "P005004", "P005006", "P005010",
                 "P015001", "P005005", "P005007", "P005008", "P005009")

chi_sf1 <- tidycensus::get_decennial(geography = "block", state = "IL",
                                     county = "Cook County", variables = vars_census,
                                     year = 2010, output = "wide", geometry = TRUE)

# Download ACS data 2011-2016
# mhhi
# % poverty
# % unemployed
# college

vars_acs <- c("C24010_004E", "C24010_040E", "C24010_001E", "B23025_005E", "B23025_003E",
              "B17021_001E", "B17021_002E", "B15002_011E", "B15002_028E", "B15002_001E",
              "B11003_016E", "B11001_001E", "B25001_001E", "B25003_003E", "B25002_002E",
              "B25002_003E", "B16004_001E", "B16004_006E", "B16004_007E", "B16004_008E",
              "B16004_011E", "B16004_012E", "B16004_013E", "B16004_016E", "B16004_017E",
              "B16004_018E", "B16004_021E", "B16004_022E", "B16004_023E", "B16004_028E",
              "B16004_029E", "B16004_030E", "B16004_033E", "B16004_034E", "B16004_035E",
              "B16004_038E", "B16004_039E", "B16004_040E", "B16004_043E", "B16004_044E",
              "B16004_045E", "B16004_050E", "B16004_051E", "B16004_052E", "B16004_055E",
              "B16004_056E", "B16004_057E", "B16004_060E", "B16004_061E", "B16004_062E",
              "B16004_065E", "B16004_066E", "B16004_067E", "B16004_008E", "B16004_013E",
              "B16004_018E", "B16004_023E", "B16004_030E", "B16004_035E", "B16004_040E",
              "B16004_045E", "B16004_052E", "B16004_057E", "B16004_062E", "B16004_067E",
              "B16004_001E", "B25038_004E", "B25038_011E", "B25038_004E", "B25038_011E",
              "B19013_001", # mhhi 
              "B15003_001", # education universe 
              "B15003_022", "B15003_023", "B15003_024", "B15003_025", # number of college educated or more
              "B25003_001", # homeowner universe
              "B25003_002") # homeowner

vars_loaded = load_variables(dataset = "acs5", year = 2013)
# View(vars_loaded)

chi_acs <- tidycensus::get_acs(geography = "block group", state = "IL",
                               county = "Cook County", variables = vars_acs,
                               year = 2013, output = "wide",
                               geometry = TRUE)

# herfindahl or fractionalization index

f_index <- function(x) 1 - rowSums(x^2)

# factor analysis

fac <- function(x, factors, rotation = "varimax", scores = "regression") {
  ok <- complete.cases(x)
  fa <- factanal(x[ok,], factors, rotation=rotation, scores=scores)
  score <- rep(NA, fa$n.obs + sum(!ok))
  score[ok] <- fa$scores
  return(score)
}

# recode Decennial Census

chi_sf1 <- chi_sf1 %>%
  transmute(
    BLK_CODE = GEOID,
    pop = P001001,
    hh = P015001,
    race_white = P005003,
    race_black = P005004,
    race_asian = P005006,
    race_hisp = P005010,
    race_other = P005005 + P005007 + P005008 + P005009,
    p_race_white = race_white / pop,
    p_race_black = race_black / pop,
    p_race_asian = race_asian / pop,
    p_race_hisp = race_hisp / pop,
    p_race_other = race_other / pop,
    age_15_35_male = P012006 + P012007 + P012008 + P012009 + P012010 +
      P012011 + P012012,
    hhi = f_index(cbind(p_race_white, p_race_black, p_race_asian,
                        p_race_hisp, p_race_other)),
    geometry
  )

chi_sf1$GEOID = substring(chi_sf1$BLK_CODE, 1, 12)

# recode ACS

chi_acs <- chi_acs %>%
  transmute(
    BG_CODE = GEOID,
    occ_management = C24010_004E + C24010_040E,
    occ_universe = C24010_001E,
    emp_unemployed = B23025_005E,
    emp_civ_lab_force = B23025_003E,
    poverty_universe = B17021_001E,
    poverty = B17021_002E,
    edu_hs = B15002_011E + B15002_028E,
    edu_universe = B15002_001E,
    hh_single_mother = B11003_016E,
    mhhi = B19013_001E, 
    hh = B11001_001E,
    hu = B25001_001E,
    hu_occupied_renter = B25003_003E,
    hu_occupied = B25002_002E,
    hu_vacant = B25002_003E,
    lang_universe = B16004_001E,
    lang_eng_limited = # Speak Eng less than "very well"
      B16004_006E + B16004_007E + B16004_008E +
      B16004_011E + B16004_012E + B16004_013E +
      B16004_016E + B16004_017E + B16004_018E +
      B16004_021E + B16004_022E + B16004_023E +
      B16004_028E + B16004_029E + B16004_030E +
      B16004_033E + B16004_034E + B16004_035E +
      B16004_038E + B16004_039E + B16004_040E +
      B16004_043E + B16004_044E + B16004_045E +
      B16004_050E + B16004_051E + B16004_052E +
      B16004_055E + B16004_056E + B16004_057E +
      B16004_060E + B16004_061E + B16004_062E +
      B16004_065E + B16004_066E + B16004_067E,
    lang_eng_no = B16004_008E + B16004_013E + B16004_018E + B16004_023E +
      B16004_030E + B16004_035E + B16004_040E + B16004_045E +
      B16004_052E + B16004_057E + B16004_062E + B16004_067E,
    lang_universe = B16004_001E,
    housing_moved_2000_2009 = B25038_004E + B25038_011E,
    housing_moved_2000_2009 = B25038_004E + B25038_011E,
    p_occ_management = occ_management / occ_universe,
    p_emp_unemployed = emp_unemployed / emp_civ_lab_force,
    p_poverty = poverty / poverty_universe,
    p_edu_hs = edu_hs / edu_universe,
    p_single_mother = hh_single_mother / hh,
    p_hu_occupied_renter = hu_occupied_renter / hu_occupied,
    p_hu_vacant = hu_vacant / hu,
    p_lang_eng_limited = lang_eng_limited / lang_universe,
    p_lang_eng_no = lang_eng_no / lang_universe,
    pcol = (B15003_022E + B15003_023E + B15003_024E + B15003_025E) / B15003_001E,
    pown = B25003_002E / B25003_001E,
    p_housing_moved_2000_2009 = housing_moved_2000_2009 / hu_occupied,
    # concentrated disadvantage, residential instability and immigrant concentration
    con_disadv = fac(cbind(p_occ_management, p_emp_unemployed, p_poverty,
                           p_edu_hs, p_single_mother), factors=1),
    res_instab = fac(cbind(p_hu_occupied_renter, p_housing_moved_2000_2009,
                           p_hu_vacant), factors=1),
    immi_con = rowMeans(cbind(p_lang_eng_limited, p_lang_eng_no))
  ) %>% dplyr::select(BG_CODE, con_disadv, res_instab, immi_con, mhhi, pcol, pown, p_poverty,
                      p_emp_unemployed)

# join decennial census with ACS data

chi_sf1$BG_CODE = substring(chi_sf1$BLK_CODE, 1, 12)
chi_bg <- chi_sf1 %>% left_join((chi_acs %>% mutate(geometry = NULL) %>% 
                                   as.data.frame), by = "BG_CODE")

# merge with boundary data

chi_shp = read_sf("geo_export_1d6332b4-f0c9-4c57-bfc8-9737febddf3d.shp") %>%
  dplyr::select(geometry)

chi_shp = st_transform(chi_shp, st_crs(chi_bg))
st_crs(chi_shp) = st_crs(chi_bg)

chi_bg = st_make_valid(chi_bg, dist = 0)
chi_shp = st_make_valid(chi_shp, dist = 0)
chi_bg_int = st_intersects(chi_shp, chi_bg)
chi_bg = chi_bg[chi_bg_int[[1]], ]

# wombling, generating difference measure

chi_bg = st_make_valid(chi_bg, dist=0)

uq_bgs = unique(chi_bg$BLK_CODE)

test = poly2nb(chi_bg)
chi_bg$p_race_black_blv = NA
chi_bg$p_race_white_blv = NA
chi_bg$p_race_hisp_blv = NA
chi_bg$p_race_asian_blv = NA

# loop to generate measure

for (i in 1:length(uq_bgs)) {
  
  print(paste0("Iteration ", i))
  
  chi_bg$p_race_black_blv[i] = 
    max(abs(chi_bg[chi_bg$BLK_CODE %in% uq_bgs[i], ]$p_race_black - 
              chi_bg[test[[i]], ]$p_race_black), na.rm = TRUE)
  chi_bg$p_race_white_blv[i] = 
    max(abs(chi_bg[chi_bg$BLK_CODE %in% uq_bgs[i], ]$p_race_white - 
              chi_bg[test[[i]], ]$p_race_white), na.rm = TRUE)
  chi_bg$p_race_hisp_blv[i] = 
    max(abs(chi_bg[chi_bg$BLK_CODE %in% uq_bgs[i], ]$p_race_hisp - 
              chi_bg[test[[i]], ]$p_race_hisp), na.rm = TRUE)
  chi_bg$p_race_asian_blv[i] = 
    max(abs(chi_bg[chi_bg$BLK_CODE %in% uq_bgs[i], ]$p_race_asian - 
              chi_bg[test[[i]], ]$p_race_asian), na.rm = TRUE)
  
}

# read in categorized 311 data
chi_311_cat <- read_csv("Chicago_311_clean.csv")

# filter to 2018-2019 calls since arrests are 2018-2023
# but first format date column
chi_311_cat<-chi_311_cat %>% separate(CREATED_DATE,sep=' ',into=c("date","time"))
chi_311_cat$date<- as.Date(chi_311_cat$date,format="%m/%d/%Y")

chi_311 = chi_311_cat %>% filter(date > "2017-12-31" & date < "2020-01-01")

# locate 311 calls within census blocks
# remove missing values from 311 call geometry
chi_311 = chi_311 %>% filter(!is.na(LATITUDE))

# convert 311 call data to sf
chi_311  <- st_as_sf(x = chi_311, coords = c("LONGITUDE", "LATITUDE"))

# set crs btw census data and 311 calls
st_crs(chi_311) = st_crs(chi_bg)
chi_311 = st_transform(chi_311, st_crs(chi_bg))

chi_311_merged <- st_join(chi_311,chi_bg,
                          join = st_intersects)

chi_311_merged = chi_311_merged %>% dplyr::select(SR_NUMBER,BLK_CODE, call_code)

# get only unique service requests
# lou_311_merged <- lou_311_merged[!duplicated(lou_311_merged$service_request_id), ]

# convert call_code to character
chi_311_merged$call_code <- as.character(chi_311_merged$call_code)

# reshape data from long form to wide form with each call code as column
chi_311_merged$val = 1
wide_311_calls = chi_311_merged %>% 
  tidyr::pivot_wider(names_from = "call_code", values_from = "val", values_fill = 0)

## sum number of each type of call per category by BG and year
wide_311_calls_sum = wide_311_calls %>% 
  group_by(BLK_CODE) %>% 
  summarise(across(c("2","8","9","11","6","10","7","4","5","12","1","13","NA","3"), sum)) %>% 
  as.data.frame %>% mutate(geometry = NULL)

# sum total calls per year
wide_311_calls_sum = wide_311_calls_sum %>%
  mutate(total_calls = rowSums(across(2:15)))

# sum top 4 categories 
wide_311_calls_sum$high_SES <- wide_311_calls_sum$`1` + wide_311_calls_sum$`7` + wide_311_calls_sum$`9` + wide_311_calls_sum$`6`
wide_311_calls_sum$high_SES_stan <- (wide_311_calls_sum$high_SES/wide_311_calls_sum$total_calls)

chi_311_final = wide_311_calls_sum

# merge with chi_bg based on BLK_CODE
chi_311_final = chi_311_final %>% dplyr::select(c(BLK_CODE,high_SES_stan))
chi_bg <- merge(chi_bg, chi_311_final, by = "BLK_CODE", all.x = TRUE)

# wombling, generating diff measure between proportion of high SES calls
chi_bg = st_make_valid(chi_bg, dist=0)

uq_bgs = unique(chi_bg$BLK_CODE)

test = poly2nb(chi_bg)
chi_bg$ses_blv = NA

# loop to generate measure

for (i in 1:length(uq_bgs)) {
  
  print(paste0("Iteration ", i))
  
  chi_bg$ses_blv[i] = 
    max(abs(chi_bg[chi_bg$BLK_CODE %in% uq_bgs[i], ]$high_SES_stan - 
              chi_bg[test[[i]], ]$high_SES_stan), na.rm = TRUE)
}

chi_bg$ses_blv = 
  ifelse(chi_bg$ses_blv == Inf | chi_bg$ses_blv == -Inf, NA, chi_bg$ses_blv)

# save chi with boundary measures
save(x = chi_bg, file = "chi_blv.RData")

### data is confidential but providing the code used to aggregate geolocated arrests by census block
### at the end of the code that is commented out, load the dataset with arrests to proceed in the script



# #### loading in chicago categorized arrest data #### 
# load("chi_arrests.RData")
# 
# 
# ## find intersection btw census data and arrest data
# library(mapview)
# library(leaflet)
# st_crs(chi_arrests) = st_crs(chi_bg)
# chi_arrests = st_transform(chi_arrests, st_crs(chi_bg))
# 
# ## find intersection btw census data and arrest data
# st_crs(chi_arrests) = st_crs(chi_bg)
# chi_arrests = st_transform(chi_arrests, st_crs(chi_bg))
# 
# chi_ints = st_intersects(chi_arrests, chi_bg)
# chi_ints2a = st_intersects(chi_arrests %>% filter(felony == 1), chi_bg)
# chi_ints2b = st_intersects(chi_arrests %>% filter(felony == 0), chi_bg)
# chi_ints2c = st_intersects(chi_arrests %>% filter(violent == 1), chi_bg)
# chi_ints2d = st_intersects(chi_arrests %>% filter(violent == 0), chi_bg)
# chi_ints2e = st_intersects(chi_arrests %>% filter(society == 1), chi_bg)
# chi_ints2f = st_intersects(chi_arrests %>% filter(person == 1), chi_bg)
# chi_ints2g = st_intersects(chi_arrests %>% filter(property == 1), chi_bg)
# 
# theout = chi_bg[chi_ints %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>% 
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(total_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2a = chi_bg[chi_ints2a %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(felony_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2b = chi_bg[chi_ints2b %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(misdemeanor_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2c = chi_bg[chi_ints2c %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(violent_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2d = chi_bg[chi_ints2d %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(nonviolent_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2e = chi_bg[chi_ints2e %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(society_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2f = chi_bg[chi_ints2f %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(person_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# theout2g = chi_bg[chi_ints2g %>% unlist, ] %>% 
#   dplyr::select(BLK_CODE) %>%
#   mutate(count = 1) %>% 
#   group_by(BLK_CODE) %>% 
#   summarize(property_arrests = sum(count, na.rm = TRUE)) %>% 
#   as.data.frame %>% 
#   mutate(geometry = NULL)
# 
# 
# chi_bg2 = merge(chi_bg, theout, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2a, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2b, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2c, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2d, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2e, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2f, by = "BLK_CODE", all.x = TRUE)
# chi_bg2 = merge(chi_bg2, theout2g, by = "BLK_CODE", all.x = TRUE)
# 
# 
# 
# chi_bg2$total_arrests = ifelse(is.na(chi_bg2$total_arrests), 0, chi_bg2$total_arrests)
# chi_bg2$felony_arrests = ifelse(is.na(chi_bg2$felony_arrests), 0, chi_bg2$felony_arrests)
# chi_bg2$misdemeanor_arrests = ifelse(is.na(chi_bg2$misdemeanor_arrests), 0, chi_bg2$misdemeanor_arrests)
# chi_bg2$violent_arrests = ifelse(is.na(chi_bg2$violent_arrests), 0, chi_bg2$violent_arrests)
# chi_bg2$nonviolent_arrests = ifelse(is.na(chi_bg2$nonviolent_arrests), 0, chi_bg2$nonviolent_arrests)
# chi_bg2$society_arrests = ifelse(is.na(chi_bg2$society_arrests), 0, chi_bg2$society_arrests)
# chi_bg2$person_arrests = ifelse(is.na(chi_bg2$person_arrests), 0, chi_bg2$person_arrests)
# chi_bg2$property_arrests = ifelse(is.na(chi_bg2$property_arrests), 0, chi_bg2$property_arrests)


### read in dataset with aggregate arrests by type to continue

load("chi_block_arrests.RData")

chi_block_arrests <- chi_bg2

 

### crime data ###
chi_crime <- read_csv("Crimes_-_2001_to_Present.csv")

# format date
chi_crime$crime_date <- chi_crime$Date
chi_crime$crime_date = substring(chi_crime$crime_date, 1, 10)
chi_crime$crime_date = as.Date(chi_crime$crime_date, format = "%m/%d/%Y")

# chi_crime$Date = as.Date(chi_crime$Date, format="%m/%d/%Y")
# chi_crime_filtered = chi_crime%>%
  # filter(Date > as.Date("2014-01-01"))%>%
  # filter(Date < as.Date("2023-01-27"))

violent_type = c("ARSON", "ASSAULT", "BATTERY", "CRIM SEXUAL ASSAULT", "CRIMINAL SEXUAL ASSAULT", 
                 "HOMICIDE", "ROBBERY", "KIDNAPPING", "RITUALISM")
violent_description = c("CRIM SEX ABUSE BY FAM MEMBER", "SEX ASSLT OF CHILD BY FAM MBR", 
                        "AGG SEX ASSLT OF CHILD FAM MBR", "AGG CRIM SEX ABUSE FAM MEMBER", 
                        "COMMERCIAL SEX ACTS", "CRIMINAL SEXUAL ABUSE", "AGG CRIMINAL SEXUAL ABUSE", 
                        "SEXUAL EXPLOITATION OF A CHILD")

property_type = c("THEFT", "BURGLARY", "CRIMINAL TRESPASS", "LIQUOR LAW VIOLATION", "ROBBERY", 
                 "MOTOR VEHICLE THEFT")

chi_crime_violent = chi_crime %>%
  mutate(violent = ifelse(`Primary Type` %in% violent_type | Description %in% violent_description, 1, 0)) %>%
  mutate(property = ifelse(`Primary Type` %in% property_type, 1, 0))

chi_crime_violent = chi_crime_violent %>% filter(!is.na(Latitude))
chi_crime_violent = chi_crime_violent%>% st_as_sf(coords = c("Longitude", "Latitude"))

st_crs(chi_crime_violent) = st_crs(chi_bg2)
chi_crime_violent = st_transform(chi_crime_violent, crs = st_crs(chi_bg2))

chi_ints3a = st_intersects(chi_crime_violent %>% filter(violent == 1), chi_bg2)
chi_ints3b = st_intersects(chi_crime_violent %>% filter(property == 1), chi_bg2)


theout3a = chi_bg2[chi_ints3a %>% unlist, ] %>% 
  dplyr::select(BLK_CODE) %>%
  mutate(count = 1) %>% 
  group_by(BLK_CODE) %>% 
  summarize(violent_crime = sum(count, na.rm = TRUE)) %>% 
  as.data.frame %>% 
  mutate(geometry = NULL)

theout3b = chi_bg2[chi_ints3b %>% unlist, ] %>% 
  dplyr::select(BLK_CODE) %>%
  mutate(count = 1) %>% 
  group_by(BLK_CODE) %>% 
  summarize(property_crime = sum(count, na.rm = TRUE)) %>% 
  as.data.frame %>% 
  mutate(geometry = NULL)


chi_bg3 = merge(chi_bg2, theout3a, by = "BLK_CODE", all.x = TRUE)
chi_bg3 = merge(chi_bg3, theout3b, by = "BLK_CODE", all.x = TRUE)

chi_bg3$crime_all = 
  chi_bg3$violent_crime + chi_bg3$property_crime


# add in com_density and phys_bound vars
load("chi_fin_pb3.RData")

# get indicator of physical boundaries by block
chi_phys = chi_fin_pb3 %>% dplyr::select(c(BLK_CODE, phys_bound))
chi_phys = chi_phys %>% as.data.frame() %>% dplyr::select(-c(geometry))

chi_bg3 <- merge(chi_bg3, chi_phys, by = "BLK_CODE")

# integrate commercial density
load("chi_cd.RData")
head(chi_cd)

# Calculate the average total jobs by census block
avg_total_jobs <- chi_cd %>%
  group_by(BLK_CODE) %>%
  summarise(total_jobs = mean(total_jobs, na.rm = TRUE))

chi_bg3 <- merge(chi_bg3, avg_total_jobs, by = "BLK_CODE")

# create com_density measure

chi_bg3$com_density <- chi_bg3$total_jobs/chi_bg3$pop
summary(chi_bg3$com_density)


#### quick edits before saving #### 

# resolving infinite numbers 

chi_bg3 = 
  chi_bg3 %>% 
  mutate(p_race_black_blv = ifelse(p_race_black_blv == Inf | 
                                     p_race_black_blv == -Inf, NA, p_race_black_blv),
         p_race_hisp_blv = ifelse(p_race_hisp_blv == Inf | 
                                    p_race_hisp_blv == -Inf, NA, p_race_hisp_blv),
         p_race_asian_blv = ifelse(p_race_asian_blv == Inf | 
                                     p_race_asian_blv == -Inf, NA, p_race_asian_blv),
         p_race_white_blv = ifelse(p_race_white_blv == Inf | 
                                     p_race_white_blv == -Inf, NA, p_race_white_blv))

# global measure

chi_bg3$edge_race_areal = 
  chi_bg3[, c("p_race_black_blv", "p_race_hisp_blv",
              "p_race_asian_blv", "p_race_white_blv")] %>% 
  mutate(geometry = NULL) %>% 
  as.data.frame %>% 
  as.matrix %>% 
  matrixStats::rowMaxs(x = .)

# pairwise measure 

chi_bg3 = 
  chi_bg3 %>% 
  mutate(edge_race_wb = p_race_black_blv * p_race_white_blv,
         edge_race_wh = p_race_hisp_blv * p_race_white_blv,
         edge_race_hb = p_race_hisp_blv * p_race_black_blv)

# arrest rate

chi_bg3$arrest_rate = (chi_bg3$total_arrests / chi_bg3$pop) * 10
chi_bg3$arrest_rate = 
  ifelse(chi_bg3$arrest_rate == Inf | chi_bg3$arrest_rate == -Inf, NA,
         chi_bg3$arrest_rate)

chi_bg3$felony_arrest_rate = (chi_bg3$felony_arrests / chi_bg3$pop) * 10
chi_bg3$felony_arrest_rate = 
  ifelse(chi_bg3$felony_arrest_rate == Inf | chi_bg3$felony_arrest_rate == -Inf, NA,
         chi_bg3$felony_arrest_rate)

chi_bg3$misdemeanor_arrest_rate = (chi_bg3$misdemeanor_arrests / chi_bg3$pop) * 10
chi_bg3$misdemeanor_arrest_rate = 
  ifelse(chi_bg3$misdemeanor_arrest_rate == Inf | chi_bg3$misdemeanor_arrest_rate == -Inf, NA,
         chi_bg3$misdemeanor_arrest_rate)

chi_bg3$violent_arrest_rate = (chi_bg3$violent_arrests / chi_bg3$pop) * 10
chi_bg3$violent_arrest_rate = 
  ifelse(chi_bg3$violent_arrest_rate == Inf | chi_bg3$violent_arrest_rate == -Inf, NA,
         chi_bg3$violent_arrest_rate)

chi_bg3$nonviolent_arrest_rate = (chi_bg3$nonviolent_arrests / chi_bg3$pop) * 10
chi_bg3$nonviolent_arrest_rate = 
  ifelse(chi_bg3$nonviolent_arrest_rate == Inf | chi_bg3$nonviolent_arrest_rate == -Inf, NA,
         chi_bg3$nonviolent_arrest_rate)

chi_bg3$society_arrest_rate = (chi_bg3$society_arrests / chi_bg3$pop) * 10
chi_bg3$society_arrest_rate = 
  ifelse(chi_bg3$society_arrest_rate == Inf | chi_bg3$society_arrest_rate == -Inf, NA,
         chi_bg3$society_arrest_rate)

chi_bg3$person_arrest_rate = (chi_bg3$person_arrests / chi_bg3$pop) * 10
chi_bg3$person_arrest_rate = 
  ifelse(chi_bg3$person_arrest_rate == Inf | chi_bg3$person_arrest_rate == -Inf, NA,
         chi_bg3$person_arrest_rate)

chi_bg3$pr